#
import unittest
import math
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, rayleigh
# from scipy.special import marcumq
from core.radar_detection import RadarDetection

class TRadarDetection(unittest.TestCase):
    # python -m unittest -v uts.core.t_radar_detection.TRadarDetection.test_fig2_2
    def test_fig2_2(self):
        # 设置参数
        mu = 0       # 高斯分布的均值
        sigma = 1.5  # 标准差

        # 创建数据点
        xg = np.linspace(-6, 6, 1500)  # 高斯分布随机变量范围
        xr = np.linspace(0, 6, 1500)   # 瑞利分布随机变量范围

        # 计算概率密度函数
        ynorm = norm.pdf(xg, mu, sigma)     # 高斯PDF
        yray = rayleigh.pdf(xr, scale=sigma) # 瑞利PDF (scale参数对应σ)

        # 绘制图形
        plt.figure(figsize=(10, 6))
        plt.plot(xg, ynorm, 'k', label='Gaussian pdf')       # 黑色实线
        plt.plot(xr, yray, 'k-.', label='Rayleigh pdf')      # 黑色点划线

        # 添加标签和图例
        plt.grid(True)
        plt.legend()
        plt.xlabel('x')
        plt.ylabel('Probability density')

        # 添加文本注释（位置根据图形手动调整）
        plt.text(-1.5, 0.25, r'$\mu = 0; \sigma = 1.5$', fontsize=12)
        plt.text(3.2, 0.15, r'$\sigma = 1.5$', fontsize=12)

        plt.tight_layout()
        plt.show()

    # python -m unittest -v uts.core.t_radar_detection.TRadarDetection.test_que_func
    def test_que_func(self):
        x = 2.2
        f = RadarDetection.que_func(x)
        print(f'value: {f};')

    # python -m unittest -v uts.core.t_radar_detection.TRadarDetection.test_fig2_3
    def test_fig2_3(self):
        # 生成虚警率对数范围 (0.01 到 250)
        logpfa = np.linspace(0.01, 250, 1000)

        # 计算变量 var = 10^(logpfa/10)
        var = 10**(logpfa / 10.0)

        # 计算归一化阈值 vtnorm = sqrt(ln(var))
        vtnorm = np.sqrt(np.log(var))

        # 创建图形
        plt.figure(figsize=(10, 6))

        # 绘制半对数图：x轴对数刻度，y轴线性刻度
        plt.semilogx(logpfa, vtnorm, 'k', linewidth=1.5)

        # 添加网格和标签
        plt.grid(True, which='both', linestyle='--', alpha=0.7)
        plt.title('Normalized Threshold vs. Log of P_{fa}')
        plt.xlabel('Log of Probability of False Alarm (log(P_{fa}))')
        plt.ylabel('Normalized Threshold (V_T / σ)')
        plt.tight_layout()

        # 显示图形
        plt.show()

    # python -m unittest -v uts.core.t_radar_detection.TRadarDetection.test_fig2_4
    def test_fig2_4(self):
        print(f'书中代码实现有问题')

    # python -m unittest -v uts.core.t_radar_detection.TRadarDetection.test_fig2_6a
    def test_fig2_6a(self):
        # Set probability values
        pfa1 = 1.0e-2
        pfa2 = 1.0e-6
        pfa3 = 1.0e-10
        pfa4 = 1.0e-13
        pd1 = 0.5
        pd2 = 0.8
        pd3 = 0.95
        pd4 = 0.999

        # Initialize arrays
        np_values = np.arange(1, 1001)
        I1 = np.zeros(len(np_values))
        I2 = np.zeros(len(np_values))
        I3 = np.zeros(len(np_values))
        I4 = np.zeros(len(np_values))

        # Calculate improvement factors
        for index, np_val in enumerate(np_values):
            I1[index] = RadarDetection.improv_fac(np_val, pfa1, pd1)
            I2[index] = RadarDetection.improv_fac(np_val, pfa2, pd2)
            I3[index] = RadarDetection.improv_fac(np_val, pfa3, pd3)
            I4[index] = RadarDetection.improv_fac(np_val, pfa4, pd4)

        # Plot results
        plt.figure(figsize=(10, 6))
        plt.semilogx(np_values, I1, 'k', label='pd=.5, nfa=e+2')
        plt.semilogx(np_values, I2, 'k--', label='pd=.8, nfa=e+6')
        plt.semilogx(np_values, I3, 'k-.', label='pd=.95, nfa=e+10')
        plt.semilogx(np_values, I4, 'k:', label='pd=.999, nfa=e+13')

        plt.xlabel('Number of pulses')
        plt.ylabel('Improvement factor I - dB')
        plt.title('Non-coherent Integration Improvement Factor')
        plt.legend()
        plt.grid(True, which="both", ls="-")
        plt.xlim(1, 1000)
        plt.ylim(0, 25)  # Adjusted based on typical output range
        plt.tight_layout()
        plt.show()

    # python -m unittest -v uts.core.t_radar_detection.TRadarDetection.test_fig2_6b
    def test_fig2_6b(self):
        # 设置参数
        pfa = 1.0e-12  # 所有情况使用相同的虚警率
        pd_values = [0.5, 0.8, 0.95, 0.99]  # 不同的检测概率

        # 生成脉冲数序列 (1-1000)
        n_P = np.arange(1, 1001)
        loss_curves = []  # 存储每条损失曲线的数据

        # 计算不同检测概率下的累积损失
        for pd in pd_values:
            L = np.zeros(len(n_P))  # 存储当前pd下的损失值
            for i, np_val in enumerate(n_P):
                # 计算改善因子 (dB)
                I = RadarDetection.improv_fac(np_val, pfa, pd)
                # 将改善因子转换为线性值
                i_linear = 10 ** (0.1 * I)
                # 计算累积损失 (dB)
                L[i] = -10 * np.log10(i_linear / np_val)
            loss_curves.append(L)

        # 设置绘图样式
        line_styles = ['k-', 'k--', 'k-.', 'k:']
        labels = [f'pd={pd}, nfa=e+12' for pd in pd_values]

        # 创建图形
        plt.figure(figsize=(10, 6))

        # 绘制所有曲线
        for i, (L, style, label) in enumerate(zip(loss_curves, line_styles, labels)):
            plt.semilogx(n_P, L, style, label=label)

        # 设置坐标轴和标签
        plt.xlabel('Number of pulses')
        plt.ylabel('Integration loss - dB')
        plt.title('Integration Loss vs. Number of Pulses')

        # 添加图例和网格
        plt.legend()
        plt.grid(True, which="both", linestyle='-', alpha=0.6)

        # 设置坐标轴范围
        plt.ylim(-1, 9)  # 根据典型值设置y轴范围
        plt.xlim(1, 1000)  # 设置x轴范围

        plt.tight_layout()
        plt.show()

    # python -m unittest -v uts.core.t_radar_detection.TRadarDetection.test_mini_case2_1
    def test_mini_case2_1(self):
        '''
        利用给定参数，完成如下任务：
        确定脉冲重复频率f_r，脉冲宽度tau，峰值功率P_t，虚警率P_fa，检出率P_D，最小可探测信号水平S_min
        采用10个非相参累积脉冲来降低传输功率需求
        采用单脉冲模式，探测距离为9公里时，估计检出率
        '''
        c = 3e8 # 光速
        f_0 = 1.5e9 # 中心频率 GHz
        B = 2e6 # 带宽 MHz
        F_dB = 8 # 噪声数 dB
        F = 10**(F_dB / 10.0)
        L_dB = 4 # 系统损失 dB
        L = 10**(L_dB / 10.0)
        T_fa = 12*60 # 虚警时间 s
        R = 12e3 # 探测距离 Km
        SNR_dB = 13.85 # 最小SNR
        SNR = 10**(SNR_dB / 10.0)
        G = 5000 # 天线增益
        sigma = 1 # RCS m^{2}
        k = 1.38e-23 # 常数
        T_e = 290.0
        # 计算脉冲重复频率
        R_u = R
        f_r = c / (2*R_u)
        print(f'脉冲重复频率：{f_r};')
        tau = 1 / B
        print(f'脉冲宽度：{tau};')
        # 计算虚警率：式2.21
        P_fa = 1 / (B * T_fa)
        print(f'虚警率为：{P_fa};')
        # 求检出率：式2.24
        P_D_alpha = np.sqrt(2*SNR)
        P_D_beta = np.sqrt(2*np.log(1/P_fa))
        print(f'alpha: {P_D_alpha}; beta: {P_D_beta};')
        P_D = RadarDetection.marcumsq_new(P_D_alpha, P_D_beta)
        print(f'检出率：{P_D};')
        lambda_ = c / f_0
        P_t = SNR * ((4*np.pi)**3 * R**4 * k * T_e * B * F * L) / (G**2 * lambda_**2 * sigma)
        print(f'峰值发射功率：{P_t} W;')
        S_min = (P_t * G**2 * lambda_**2 * sigma) / ((4*np.pi)**3 * R**4 * L)
        print(f'最小信号（电压）：{S_min} V;')
        # 2.
        n_P = 10
        I_n_P = RadarDetection.improv_fac(n_P, P_fa, P_D)
        I_n_P_dB = 10*np.log10(I_n_P)
        print(f'非相参累积改善因子：{I_n_P}; or {I_n_P_dB}dB;')
        L_NCI = n_P / I_n_P
        L_NCI_dB = 10**np.log10(L_NCI)
        print(f'累积损失：{L_NCI}; {L_NCI_dB}dB;')
        SNR_NCI_dB = SNR_dB - I_n_P_dB + L_NCI_dB
        print(f'SNR_NCI: {SNR_NCI_dB};')

    # python -m unittest -v uts.core.t_radar_detection.TRadarDetection.test_fig2_7
    def test_fig2_7(self):
        # Main program to reproduce the plot
        x = np.arange(0, 20.1, 0.1)  # From 0 to 20 in steps of 0.1
        val1 = []
        val2 = []
        val3 = []
        val4 = []

        for xi in x:
            val1.append(RadarDetection.incomplete_gamma(xi, 1))
            val2.append(RadarDetection.incomplete_gamma(xi, 3))
            val3.append(RadarDetection.incomplete_gamma(xi, 6))
            val4.append(RadarDetection.incomplete_gamma(xi, 10))

        plt.plot(x, val1, 'k', label='N = 1')
        plt.plot(x, val2, 'k:', label='N = 3')
        plt.plot(x, val3, 'k--', label='N = 6')
        plt.plot(x, val4, 'k-.', label='N = 10')
        plt.legend()
        plt.xlabel('x')
        plt.ylabel('Incomplete Gamma function (x,N)')
        plt.grid(True)
        plt.show()

    # python -m unittest -v uts.core.t_radar_detection.TRadarDetection.test_fig2_8
    def test_fig2_8(self):
        n_values = list(range(1, 151))  # 脉冲数从1到150
        y1, y3, y4 = [], [], []  # 存储不同nfa的阈值
        
        for n in n_values:
            _, vt1 = RadarDetection.threshold(1000, n)
            y1.append(vt1)
            
            _, vt2 = RadarDetection.threshold(10000, n)
            y3.append(vt2)
            
            _, vt3 = RadarDetection.threshold(500000, n)
            y4.append(vt3)
        
        # 创建对数坐标图
        plt.figure(figsize=(10, 6))
        plt.loglog(n_values, y1, 'k', label='nfa=1000')
        plt.loglog(n_values, y3, 'k--', label='nfa=10000')
        plt.loglog(n_values, y4, 'k-.', label='nfa=500000')
        
        plt.axis([1, 200, 1, 300])  # 设置坐标轴范围
        plt.xlabel('Number of pulses')
        plt.ylabel('Threshold')
        plt.legend()
        plt.grid(True, which="both", ls="-")
        plt.title('Threshold vs Number of Pulses')
        plt.show()

    # python -m unittest -v uts.core.t_radar_detection.TRadarDetection.test_fig2_9
    def test_fig2_9(self):
        pfa = 1e-9
        nfa = math.log(2) / pfa
        b = math.sqrt(-2.0 * math.log(pfa))
        
        snr_db = np.arange(0, 20.1, 0.1)  # SNR从0到20dB，步长0.1dB
        pro = []       # 单脉冲检测概率
        prob205 = []   # 10个脉冲累积的检测概率
        
        # 计算不同SNR下的检测概率
        for snr in snr_db:
            a = math.sqrt(2.0 * 10**(0.1 * snr))
            pro.append(RadarDetection.marcumsq_new(a, b))
            prob205.append(RadarDetection.pd_swerling5(pfa, 1, 10, snr))
        
        # 绘制结果
        plt.figure(figsize=(10, 6))
        plt.plot(snr_db, pro, 'k', label='np = 1')
        plt.plot(snr_db, prob205, 'k:', label='np = 10')
        
        plt.axis([0, 20, 0, 1])
        plt.xlabel('SNR - dB')
        plt.ylabel('Probability of detection')
        plt.title('Detection Probability for Swerling V Targets')
        plt.legend()
        plt.grid(True)
        plt.show()

    # python -m unittest -v uts.core.t_radar_detection.TRadarDetection.test_fig2_10
    def test_fig2_10(self):
        # Reproduce Fig. 2.10
        pfa = 1e-9
        nfa = np.log(2) / pfa
        b = np.sqrt(-2.0 * np.log(pfa))

        snr_values = np.arange(0, 22.1, 0.1)
        pro = []  # Swerling V (non-fluctuating)
        prob = [] # Swerling I

        index = 0
        for snr in snr_values:
            index += 1
            a = np.sqrt(2.0 * 10**(0.1 * snr))
            pro.append(RadarDetection.marcumsq_new(a, b))  # Marcum Q-function for Swerling V
            prob.append(RadarDetection.pd_swerling1(nfa, 1, snr))

        # Plot results
        plt.figure()
        plt.plot(snr_values, pro, 'k', label='Swerling V')
        plt.plot(snr_values, prob, 'k:', label='Swerling I')
        plt.axis([2, 22, 0, 1])
        plt.xlabel('SNR (dB)')
        plt.ylabel('Probability of Detection')
        plt.legend()
        plt.grid(True)
        plt.show()

    # python -m unittest -v uts.core.t_radar_detection.TRadarDetection.test_fig2_11ab
    def test_fig2_11ab(self):
        # 主程序 - 重现图2.11
        pfa = 1e-8
        self._fig2_11(pfa=pfa)
        pfa = 1e-11
        self._fig2_11(pfa=pfa)

    def _fig2_11(self, pfa):
        nfa = math.log(2) / pfa
        # 创建SNR范围
        snr_values = np.arange(-10, 30.5, 0.5)
        # 初始化概率数组
        prob1 = []    # np = 1
        prob10 = []   # np = 10
        prob50 = []   # np = 50
        prob100 = []  # np = 100
        # 计算检测概率
        for snr in snr_values:
            prob1.append(RadarDetection.pd_swerling1(nfa, 1, snr))
            prob10.append(RadarDetection.pd_swerling1(nfa, 10, snr))
            prob50.append(RadarDetection.pd_swerling1(nfa, 50, snr))
            prob100.append(RadarDetection.pd_swerling1(nfa, 100, snr))

        # 绘图
        plt.figure(figsize=(12, 8))

        plt.plot(snr_values, prob1, 'k-', linewidth=2, label='np = 1')
        plt.plot(snr_values, prob10, 'k:', linewidth=2, label='np = 10')
        plt.plot(snr_values, prob50, 'k--', linewidth=2, label='np = 50')
        plt.plot(snr_values, prob100, 'k-.', linewidth=2, label='np = 100')

        plt.axis([-10, 30, 0, 1])
        plt.xlabel('SNR (dB)', fontsize=12)
        plt.ylabel('Probability of Detection', fontsize=12)
        plt.title('Swerling I Detection Probability for Different Pulse Numbers', fontsize=14)
        plt.legend(loc='lower right', fontsize=10)
        plt.grid(True, linestyle='--', alpha=0.7)

        # 添加子图标识
        plt.text(0.05, 0.95, 'Fig. 2.11', transform=plt.gca().transAxes, 
                fontsize=12, verticalalignment='top')

        plt.tight_layout()
        plt.show()
        
    # python -m unittest -v uts.core.t_radar_detection.TRadarDetection.test_fig2_12
    def test_fig2_12(self):
        pfa = 1e-10
        nfa = math.log(2) / pfa
        # 创建SNR范围
        snr_values = np.arange(-10, 30.5, 0.5)
        # 初始化概率数组
        prob1 = []    # np = 1
        prob10 = []   # np = 10
        prob50 = []   # np = 50
        prob100 = []  # np = 100
        # 计算检测概率
        for snr in snr_values:
            prob1.append(RadarDetection.pd_swerling2(nfa, 1, snr))
            prob10.append(RadarDetection.pd_swerling2(nfa, 10, snr))
            prob50.append(RadarDetection.pd_swerling2(nfa, 50, snr))
            prob100.append(RadarDetection.pd_swerling2(nfa, 100, snr))

        # 绘图
        plt.figure(figsize=(12, 8))

        plt.plot(snr_values, prob1, 'k-', linewidth=2, label='np = 1')
        plt.plot(snr_values, prob10, 'k:', linewidth=2, label='np = 10')
        plt.plot(snr_values, prob50, 'k--', linewidth=2, label='np = 50')
        plt.plot(snr_values, prob100, 'k-.', linewidth=2, label='np = 100')

        plt.axis([-10, 30, 0, 1])
        plt.xlabel('SNR (dB)', fontsize=12)
        plt.ylabel('Probability of Detection', fontsize=12)
        plt.title('Swerling I Detection Probability for Different Pulse Numbers', fontsize=14)
        plt.legend(loc='lower right', fontsize=10)
        plt.grid(True, linestyle='--', alpha=0.7)

        # 添加子图标识
        plt.text(0.05, 0.95, 'Fig. 2.11', transform=plt.gca().transAxes, 
                fontsize=12, verticalalignment='top')

        plt.tight_layout()
        plt.show()

        
        
    # python -m unittest -v uts.core.t_radar_detection.TRadarDetection.test_fig2_13
    def test_fig2_13(self):
        pfa = 1e-9
        nfa = math.log(2) / pfa
        # 创建SNR范围
        snr_values = np.arange(-10, 30.5, 0.5)
        # 初始化概率数组
        prob1 = []    # np = 1
        prob10 = []   # np = 10
        prob50 = []   # np = 50
        prob100 = []  # np = 100
        # 计算检测概率
        for snr in snr_values:
            prob1.append(RadarDetection.pd_swerling3(nfa, 1, snr))
            prob10.append(RadarDetection.pd_swerling3(nfa, 10, snr))
            prob50.append(RadarDetection.pd_swerling3(nfa, 50, snr))
            prob100.append(RadarDetection.pd_swerling3(nfa, 100, snr))

        # 绘图
        plt.figure(figsize=(12, 8))

        plt.plot(snr_values, prob1, 'k-', linewidth=2, label='np = 1')
        plt.plot(snr_values, prob10, 'k:', linewidth=2, label='np = 10')
        plt.plot(snr_values, prob50, 'k--', linewidth=2, label='np = 50')
        plt.plot(snr_values, prob100, 'k-.', linewidth=2, label='np = 100')

        plt.axis([-10, 30, 0, 1])
        plt.xlabel('SNR (dB)', fontsize=12)
        plt.ylabel('Probability of Detection', fontsize=12)
        plt.title('Swerling I Detection Probability for Different Pulse Numbers', fontsize=14)
        plt.legend(loc='lower right', fontsize=10)
        plt.grid(True, linestyle='--', alpha=0.7)

        # 添加子图标识
        plt.text(0.05, 0.95, 'Fig. 2.11', transform=plt.gca().transAxes, 
                fontsize=12, verticalalignment='top')

        plt.tight_layout()
        plt.show()

        
        
    # python -m unittest -v uts.core.t_radar_detection.TRadarDetection.test_fig2_14
    def test_fig2_14(self):
        pfa = 1e-9
        nfa = math.log(2) / pfa
        # 创建SNR范围
        snr_values = np.arange(-10, 30.5, 0.5)
        # 初始化概率数组
        prob1 = []    # np = 1
        prob10 = []   # np = 10
        prob50 = []   # np = 50
        prob100 = []  # np = 100
        # 计算检测概率
        for snr in snr_values:
            prob1.append(RadarDetection.pd_swerling4(nfa, 1, snr))
            prob10.append(RadarDetection.pd_swerling4(nfa, 10, snr))
            prob50.append(RadarDetection.pd_swerling4(nfa, 50, snr))
            prob100.append(RadarDetection.pd_swerling4(nfa, 100, snr))

        # 绘图
        plt.figure(figsize=(12, 8))

        plt.plot(snr_values, prob1, 'k-', linewidth=2, label='np = 1')
        plt.plot(snr_values, prob10, 'k:', linewidth=2, label='np = 10')
        plt.plot(snr_values, prob50, 'k--', linewidth=2, label='np = 50')
        plt.plot(snr_values, prob100, 'k-.', linewidth=2, label='np = 100')

        plt.axis([-10, 30, 0, 1])
        plt.xlabel('SNR (dB)', fontsize=12)
        plt.ylabel('Probability of Detection', fontsize=12)
        plt.title('Swerling I Detection Probability for Different Pulse Numbers', fontsize=14)
        plt.legend(loc='lower right', fontsize=10)
        plt.grid(True, linestyle='--', alpha=0.7)

        # 添加子图标识
        plt.text(0.05, 0.95, 'Fig. 2.11', transform=plt.gca().transAxes, 
                fontsize=12, verticalalignment='top')

        plt.tight_layout()
        plt.show()

    # python -m unittest -v uts.core.t_radar_detection.TRadarDetection.test_fluct_loss
    def test_fluct_loss(self):
        P_D, pfa, n_P, swerling_type = 0.65, 1e-9, 10, 1
        Lf, PD_SW5, SNR_SW5, PD_SW, SNR_SW = RadarDetection.fluct_loss(P_D, pfa, n_P, swerling_type)
        print(f'第5类：{PD_SW5}; {SNR_SW5};')
        print(f'第1类：{PD_SW}; {SNR_SW};')
        print(f'Lf: {Lf};')

    # python -m unittest -v uts.core.t_radar_detection.TRadarDetection.test_fig2_15
    def test_fig2_15(self):
        # 初始化空列表存储结果
        Lf1 = []
        Lf4 = []
        pd_values = np.arange(0.01, 1.01, 0.05)  # 从0.01到1（包含1），步长0.05

        # 遍历每个检测概率值
        for pd in pd_values:
            # 计算 Swerling I & II 的波动损失
            Lf, _, _, _, _ = RadarDetection.fluct_loss(pd, 1e-9, 1, 1)
            Lf1.append(Lf)
            
            # 计算 Swerling III & IV 的波动损失
            Lf, _, _, _, _ = RadarDetection.fluct_loss(pd, 1e-9, 1, 4)
            Lf4.append(Lf)

        # 绘制图形
        plt.figure()
        plt.plot(pd_values, Lf1, 'k', label='Swerling I & II')        # 黑色实线
        plt.plot(pd_values, Lf4, 'k:', label='Swerling III & IV')     # 黑色虚线
        plt.xlabel('Probability of detection')
        plt.ylabel('Fluctuation loss - dB')
        plt.legend()
        plt.title('Pfa=1e-9, np=1')
        plt.grid(True)
        plt.show()

    # python -m unittest -v uts.core.t_radar_detection.TRadarDetection.test_mini2_2
    def test_mini2_2(self):
        P_D = 0.5
        P_fa = 1e-7
        # 求出信噪比
        snr = -math.log(P_fa)
        snr_dB = 10*np.log10(snr)
        print(f'snr: {snr}; {snr_dB}dB;')
        R = 16e3
        SNR_R = snr_dB - 40*math.log10(R) + 40*math.log10(10e3)
        print(f'SNR_R: {SNR_R};')

